/*
 * fft_radix4.hpp
 *
 *  Created on: May 20, 2010
 *      Author: z1083248
 */

#if !defined(FFT_RADIX4_HPP_)
#define FFT_RADIX4_HPP_

#include "fft_constants.hpp"
#include "fft_ops.hpp"
#include <complex>

namespace FftLib
{
namespace BasicRadix2
{

using namespace FftLib::Common;

class Fft
{

public:
   //
   // Decimation in time (DIT) radix-2 FFT
   // Depth-first version:
   // -- optimal memory access
   // -- many more trig computations (can be placed in lookup table)
   //
   // ldn = log2(n) == 12 (for n = 4096 [64x64]);
   //
   void decimationInTimeRadix2DepthFirst4096(std::complex<float> *f,
         uint32 ldn, int is)
   {
      const uint32 n = 1UL << ldn; // 4096
      const float32 pi = is * M_PI;

      for (uint32 ldm = 1; ldm <= ldn; ++ldm) // [1, 2, 3, ..., 12]
      {
         const uint32 m = (1UL << ldm); // [2, 4, 8, 16, ..., 4096]
         const uint32 mh = (m >> 1); // [1, 2, 4, 8, ..., 2048]
         const float32 phi = pi / (mh);

         for (uint32 r = 0; r < n; r += m) // 0 < 4096 step [2, 4, 8, 16, ..., 4096]
         {
            for (uint32 j = 0; j < mh; ++j) // 0 < [1, 2, 4, 8, ..., 2048]
            {
               uint32 i0 = r + j;
               uint32 i1 = i0 + mh;

               std::complex<float32> u = f[i0];
               std::complex<float32> v = f[i1] * SinCos(phi * j);

               f[i0] = u + v;
               f[i1] = u - v;
            }
         }
      }
   }

   void decimationInTimeRadix2DepthFirst2d4096(std::complex<float> *f,
         uint32 ldn, int is)
   {
      //      // do rows
      //      decimationInTimeRadix2DepthFirst4096(f, ldn, is);
      //      // do columns
      //      squareMatrixTransposeComplex(f, const4096::WindowWidth);
      //      decimationInTimeRadix2DepthFirst4096(f, ldn, is);
      //      squareMatrixTransposeComplex(f, const4096::WindowWidth);

      // Rows
      for (uint32 index = 0; index < const4096::WindowArea; index
            += const4096::WindowWidth)
      {
         decimationInTimeRadix2DepthFirst4096(f + index, 6, is);
      }

      // Columns
      squareMatrixTransposeComplex(f, const4096::WindowWidth);
      for (uint32 index = 0; index < const4096::WindowArea; index
            += const4096::WindowWidth)
      {
         decimationInTimeRadix2DepthFirst4096(f + index, 6, is);
      }
      squareMatrixTransposeComplex(f, const4096::WindowWidth);
   }

};

}
}

#endif /* FFT_RADIX4_HPP_ */
